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This  papier  investigates  certain  properties 
of  a method  which  calculates  a sample  of  a 
stationary  Gaussian  random  process  with  a 
specified  power  spectral  density  function. 
The  study  determines  to  what  degree  this 
method  simulates  a linear  filter.  Also 


included  are  correlation  analyses  of  equi- 
distributed  sequences  which  are  used  in  the 
method. 


I 

I 

l 

1 

I 


1 


T 

1 


TABLE  OF  CONTENTS 


Chapter 

ABSTRACT 

I.  INTRODUCTION 

II.  FRANKLIN'S  METHOD  - A FILTER 

III.  COMPUTER  PROGRAMS 

IV.  EXAMPLE  RUNS  AND  RESULTS 
APPENDIX 

BIBLIOGRAPHY 


Page 

ii 

1 

k 

10 

20 

28 

36 


I.  INTRODUCTION 


For  r a random  variable,  a collection  of  functions  of  time 
x(t)  = xr(t)  is  called  a random  process.  If  for  every  finite  collection  of 
times  t^  < • • * < t^,  the  random  variables  x(t^),  x(t2),  • • • , x('tn)  have 
a multivariate  Gaussian  distribution,  the  process  is  called  Gaussian.  The 
process  is  called  stationary  if,  for  any  increment  At,  the  random  variables 
x(t^  + At)  and  x(t^)  have  the  same  joint  distribution.  A stationary 
Gaussian  random  process  is  called  a time  series  if  t represents  values  of 
time.  1 

Mr.  Joel  N.  Franklin  proposed  a method  for  computing  a time  series 

2 

with  given  statistical  properties.  Franklin's  method  constructs  a sample 
of  a stationary  Gaussian  random  process  with  a given  power  spectral  density 
function.  This  technique  requires  the  use  of  a sequence  of  random  numbers. 
Walter  Matuska  wrote  a computer  program  for  the  Control  Data  Corporation 
3200  computer  which  produces  these  random  numbers.^  Matuska  also  wrote  a 
computer  program  GAUSSIAN  using  the  output  of  the  random  number  generator 
which  applies  Franklin's  technique  in  computing  a time  series  with  a 


■'"J.  N.  Franklin,  "Numerical  Simulation  of  Stationary  and  Nonstationary 
Gaussian  Random  Processes,"  SIAM  Review  J_,  No.  1,  p.  68  (January,  1965). 


2Ibid. , pp.  68  - 80. 

A.  Matuska,  Jr.,  and  G.  S.  Innis,  "Generation  of  a Stationary  Gaussian 
Random  Process  with  a Specified  Power  Spectral  Density  Function,”  Defense 
Research  Laboratory,  The  University  of  Texas,  Acoustical  Report  A-258 
(6  July  1966),  pp.  12  - 16. 
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specified  spectrum.  Still  another  program  SPECu  uses  this  time  series  as 
input  and  computes  the  power  spectral  density  function. 

The  applications  of  Franklin's  theory  are  varied.  One  problem  is 
that  of  detecting  a known  signal  in  various  noise  backgrounds.  The  method 
presented  in  this  paper  generates  a time  series  which  has  the  same  power 
spectrum  as  the  time  series  obtained  from  a specified  noise  background.  To 
this  time  series  can  be  added  the  time  series  for  a signal  with  known  power 
spectrum.  By  varying  the  power  spectral  density  function  of  the  backgroi*nd, 
the  effect  of  different  noise  backgrounds  upon  the  detection  technique  ^an 
be  studied.^  This  method  can  also  be  used  in  the  theory  of  optimum  linear 
prediction  and  filtering.  The  theory  concerns  the  designing  of  engirjeering 
systems  that  either  project  into  the  future  by  information  obtained  in  the 
past  or  recover  desired  signals  which  have  been  distorted  by  random  noise. 
These  systems  can  be  applied  to  communications,  meteorological  forecasting, 
and  economic  analysis.  An  example  is  the  problem  of  designing  an  optimum 
filter  to  smooth  a perturbed  message  (signal  plus  noise).  Separate  and 
independent  information  about  a second  perturbed  message  whose  signal  portion 
is  related  to  the  first  is  then  employed.  Franklin's  method  can  be  used  to 
obtain  a time  series  with  known  power  spectral  density  function  for  both 
perturbed  signals.  Then  using  known  signal  relationships,  the  time  series 

If 

Ibid. , pp.  20  - 29. 

<5 

G.  E.  Ellis  and  R.  L.  Boston,  "Program  SPECT:  Power  Spectral  Analysis  of 
a Random  Process  Using  the  CDC  1604  Computer,"  Defense  Research  Laboratory, 
The  University  of  Texas,  Austin,  Texas,  Acoustical  Report  A-252 
(December  1965)* 

^Matuska,  o£.  clt. , p.  3> 
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of  the  noise  can  be  obtained.  One  further  application  is  the  study  of 
smoothing  techniques  such  as  hamming  and  hanning  windows  in  which  it  is  of 
interest  to  know  how  the  minor  lobes  of  a power  spectral  density  function 
are  affected  by  these  smoothing  techniques.  This  can  be  accomplished  by 
constructing  time  series  for  which  the  distances  between  the  major  and 
minor  lobes  of  their  power  spectra  vary.  The  relative  sizes  of  these  lobes 
can  also  be  varied;  therefore,  the  effect  of  a smoothing  technique  on  a 

g 

given  power  spectral  density  function  can  be  studied. 

This  paper  uses  the  outputs  of  the  aforementioned  computer  programs 
to  determine  the  degree  to  which  Franklin's  method  is  a linear  filter.  The 
task  is  partially  accomplished  with  the  use  of  a computer  program  FOURXFRM 
written  for  the  Control  Data  3200  computer.  This  program  computes  the 
Fourier  transform  of  a given  time  series.  The  paper  also  includes  corre- 
lation analyses  which  are  performed  to  determine  the  degree  to  which  the 
random  number  generator  produces  independent  sequences. 


7 

J.  S.  Bendat,  Principles  and  Applications  of  Random  Noise  Theory 
(John  Wiley  and  Sons,  New  York,  195^),  pp.  201  - 206. 

Q 

R.  B.  Blackman  and  J.  W.  Tukey,  The  Measurement  of  Power  Spectra 
(Dover  Publications,  Inc.,  New  York,  1958) , PP-  33  - 37 • 
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II.  FRANKLIN'S  METHOD  - A FILTER 


To  compute  a sample  of  a stationary  Gaussian  random  process  by 
Franklin's  Method,  the  autocorrelation  function  or  the  power  spectral 
density  function  must  be  known.  The  autocorrelation  function  is  defined  as 


R(t1,t2)  = E x(tL)  x(t2) 


where  E 


x. 

1 


, the  expected  value  operator,  is  defined  as 

N 


= lim 

N-xx) 


I Ixi 


i=l 


o 

if  the  limit  exists.  If  x(t)  is  a time  series,  then  the  autocorrelation 
function  can  be  written  as 


R(t)  = E 


x(t ) x(t--r) 


If  V(oj)  denotes  the  power  spectral  density  function  for  a given 
frequency  co,  then  V(oj)  is  defined  by  first  defining 


C(cv)  = 


rT  /.  \ -icut 
/ x(t)e  dt 

J_rj] 


for  a time  series  x(t),  |t|  ^ T.  Then  let 


/Vn 

-®C  +“ 


P (cc  , X,  ten)  = / c(co)do)  , 


c 2 


where  P^o^,  x,  £*o)  measures  the  average  power  for  a range  of  frequencies 


^Franklin,  o|>.  cit. , p.  68 


h 


>m»ii  i h»ihiiwt  ‘mu  <1  * mm 
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of  width  £tx>  and  centered  at  go  . Then  V(go)  can  be  defined 

c 


V(go)  = lim 

&o-»o 


lim  P (go  , x,  Zko) 
T->“> 


If  x(t)  is  a stationary  random  process,  V(go)  can  be  shown  to  be 


V(go)  = / R(f)e  iQ)T  dr 


(2.1) 


if  the  integral  exists.  That  is,  the  power  spectral  density  function  is 
the  Fourier  transform  of  the  autocorrelation  function. 

If  x(t)  is  the  input  time  series  to  a linear  filter  and  y(t)  is 
the  output  time  series,  then  they  are  related  by  the  convolution 


U 

y(t)  = J g(t-s)  x(s)  ds 


(2.2) 


If  x is  equal  to  the  Dirac  delta  function,  then  y(t)  is  the  response  of  the 

linear  filter  to  a delta  function  input  and  equals  g(t).  It  can  then  be 

shown  that  the  power  spectral  density  functions  V (go)  and  V (go)  of  the  out- 

y x 

put  and  input  are  related  by  the  equation 


V (o>)  = |G(co)|  Vx(go)  , 


(2.3) 


where  G(oo)  is  the  Fourier  transform  of  g(t).' 


Bendat,  o£.  cit . , pp.  44,  66  - 67. 

L1W.  B.  Davenport  and  W.  L.  Root,  An  Introduction  to  the  Theory  of  Random 
Signals  and  Noise  (McGraw-Hill,  New  York,  195^),  pp.  182  - I03,  225  - 227- 


The  power  spectral  density  Vx(cd)  of  a random  process  x(t)  is 


defined  to  be  one  for  all  real  oo.  Such  a random  process  is  known  as  white 


noise.  If  such  a process  is  obtained  then  from  Eq.  (2.3), 


V(oo)  = V (a>)  = |G(a>)|‘ 


(2.4) 


It  has  been  shown  by  Davenport  and  Root  that  if 


0 < V(u))  <»,  V(o>)  = V(-to),  and  V(to)  -»  0 as  co  -»  ±°°,  (2.5) 


and  if  V(cn)  can  be  represented  as  a quotient  of  two  polynomials  in  cn,  then 


»<->-  SHSf 


(2.6) 


for  a)  real.  P and  Q are  polynomials  with  real  coefficients,  where  the 


degree  of  P is  less  than  the  degree  of  Q and  where  the  zeros  of  Q(z)  lie  in 


the  half  plane  Re(z)  < 0.  Now  G(cd)  can  be  defined 


°<»>  - 


(2.7) 


if  D denotes  the  differential  operator,  — , 


y(t)  ■ 


(2.8) 


The  next  step  in  obtaining  the  specified  time  series  is  to  solve  the 


differential  equation 


Q(D)  $(t)  = x(t)  , -oo  < t < oo 


(2.9) 


Once  Eq.  (2.9)  has  been  solved,  the  solution  of  4>(t)  is  used  to  obtain  y(t) 


by  combining  the  derivatives  of  ®(t)  of  order  lower  than  the  degree  of  Q 


y(t)  = P(D)  ®(t) 


(2.10) 


Tbid.,  pp.  232  - 234,  375  - 376. 


To  obtain  a solution  to  Eq.  (2.10),  a white  sequence  must  first  be 
generated. ^ If  the  white  sequence  is  denoted  x^,  a sequence  w^  (called  the 

w-sequence)  of  independent  samples  from  a Gaussian  distribution  is  constructed 

l. 

W2n-1  = ln  X2n-l^2cOS  2rtX2n  » and 


W2n  = ("2  ln  x2n-l)2sin  2ltx2n 


for  n = 1,  2,  3, 


When  this  w-sequence  is  used  in  the  solution  of 


(2.9))  n- dimensional  vectors  are  formed  from  successive  elements  of  the 


sequence 


where  n is  the  degree  of  Q.  As  can  be  seen  from  this  definition,  n elements 
of  the  w-sequence  are  required  for  each  vector  defined.  For  this  reason  the 
shaped  time  series  has  only  times  the  number  of  points  as  the  input  time 
series.  Here  it  is  of  interest  to  consider  the  unpublished  work  of 
Mr.  Joe  England.  For  the  same  w-sequence  England  defined  the  vectors 


^Matuska,  Op.  cit . , pp.  12  - 1 6. 
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Using  this  definition  the  shaped  time  series  is  reduced  hy  only  (n-l)  points, 
and  there  seems  to  be  no  weakening  of  the  power  of  the  method  when  these 
vectors  are  used. 

Numerically  there  does  not  exist  a white  noise  sequence  with  power 
spectral  density  exactly  one.  A correlation  analysis  of  these  "white" 
sequences  is  presented  in  the  Appendix. 

The  next  step  is  that  of  investigating  properties  of  Franklin's 
method  as  a linear  filter.  For  any  linear  filter  let  x(t)  be  the  input 
time  series  and  y(t)  the  output  time  series  such  that  they  are  related  as 
in  Eq.  (2.2). 

Theorem  I 

If  x(t)  = 0 for  t * 0 and 


:(t)dt  = 1 , 


then  y (t)  = h(t),  the  impulse  response.  If  x(t)  is  any  input,  then 


p 00 

yx(t)  = I x(t)  h(x-t)  dt 


Theorem  II 


If  V (co)  is  the  power  spectral  density  function  of  y (t)  and  V (cd) 
V x X 

^ X 

is  the  power  spectral  density  function  of  x(t),  then 

V (cd)  = |h(cd)|2V  (cd)  , 
yx 

where  H(od)  is  equal  to  the  Fourier  transform  of  h(t). 

ist  be 
1 P(io> 


The  following  question  now  must  be  answered.  Does 

2 


I H(co)  | 


Q(io>) 


when  considered  in  Franklin's  terminology? 


The  first  step  toward  obtaining  the  answer  to  the  question  is  to 
measure  the  impulse  response  of  a particular  filter.  This  is  accomplished 


by  use  of  Matuska's  program  GAUSSIAN.  From  the  program  for  an  input  time 
sequence  with  a Gaussian  distribution  and  for  a certain  ratio  of  polynomials, 
a resultant  time  series  is  obtained.  For  the  purposes  of  this  paper  the  in- 
put sequence  was  replaced  by  an  impulse  sequence  defined  to  be  zero  for  all  t 
except  for  one  point  at  which  it  is  defined  to  be  one.  With  this  sequence  as 
input,  the  resultant  time  series  approximates  the  impulse  response. 

From  the  two  theorems  presented  it  can  be  seen  that  if 

|H(»)|2  - IlMf  < 

then  | H(co) | must  equal  the  power  spectral  density  function  for  the  filter. 

15 

The  power  spectrum  is  obtained  from  the  program  SPECT.  This  equivalence 
is  shown  using 

V (cd)  = |H(a>)|2V(a>) 
yx  X 

obtained  from  Theorem  II  and  Eq.  (2.6).  For  a white  sequence  Vx(co)  is 
defined  as  onej  therefore 


Vy  <“>  ' 

X 


III.  COMPUTER  PROGRAMS 


This  discussion  concerns  three  programs  written  for  the  Control 
Data  3200  computer.  GAUSSIAN  was  written  by  W.  Matuska  to  shape  a 

16 

w-sequence  into  a time  series  with  a given  power  spectral  density  function. 

The  input  w-sequence  is  a sequence  of  independent  samples  from  a Gaussian 

distribution.  The  Gaussian  distribution  is  defined  from  a "white  sequence;" 

that  is,  an  equidistributed  sequence  7^  on  the  range  0 < 7^  < 1.  The  given 

power  spectral  density  function  V(co)  of  the  output  time  series  must  be  a 

rational  function  satisfying  the  conditions  (2.5).  One  other  input  to 

GAUSSIAN  is  the  time  increment  At  at  which  samples  of  the  specified  time 

series  are  to  be  calculated.  The  resultant  time  series  is  buffered  out  on 

17 

tape  to  be  used  as  input  to  the  program  SPECT. 

The  program  SPECT  computes  the  autocorrelation  function  of  the 
time  series  shaped  by  GAUSSIAN.  Then  the  power  spectral  density  function 
(which  is  the  Fourier  transform  of  the  autocorrelation  function)  is 
computed.  The  autocorrelation  function  is  plotted  linearly  against  time, 
and  the  power  spectral  density  function  is  plotted  against  frequency  on  a 
log- log  plot. 

In  addition  to  use  as  it  was  written,  GAUSSIAN  was  also  modified 
to  obtain  the  impulse  response  for  a particular  ratio  of  polynomials.  This 
was  accomplished  simply  by  redefining  the  input  w-sequence  to  be  an  impulse 
sequence.  The  impulse  sequence  is  defined  to  be  zero  everywhere  except  for 
one  t. 

1 Matuska,  oj>.  eft.,  pp.  20  - 29. 

^Ellis,  o£.  eft. 
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Another  program  FOURXFRM  was  written  to  compute  the  Fourier 
transform  of  the  impulse  response  as  obtained  from  GAUSSIAN.  The  program 
buffers  the  data  from  a magnetic  tape  and  computes  the  Fourier  transform 
for  a given  increment  of  time  At  and  for  a specified  frequency  range  and 
frequency  increment. 

Letting  f(t)  be  a time  series,  it  can  be  represented  as  a Fourier 
- 18 


series  as  follows: 


f(t)  =(r  / [AM  cos(cot)  + BM  sin  (cot)]  duo  > 


where 


A M = / cos  (cot ) f(t)dt  , and 


BM  = / sin(cot)  f(t)dt 


Let  t^  and  t^  be  two  consecutive  points  in  time;  then  the 
contribution  to  A(co)  from  this  interval  is 

r f(t  ) - f(t  ) 

AA(co)  = / cos  (cot ) f(t^)  + (t-t^)  dt 

4-  *■  2 1 


f(tg)  - f(tx)  p2 


- At  ; p 2 f(tjt  - t f(t  ) rx2 

—  / t cos  (cot)  dt  + 1 / cos  (cot)  dt 

" ti  J ' tl  J 


f(t2)  - f(t  ) r 1 f(t2)sin(ojt2)  - f (t^sinM^ 

= cos  (cot  ) - cos  (cot  ) + 

<t2  - tl>»2  2 1 J 


H.  Herbert  Howe,  "Fourier  Analysis  of  Nonperiodic  Pulses  on  Automatic 
Computers,"  U.  S.  Department  of  Commerce,  National  Bureau  of  Standards 
Report  5018  (5  October  1956),  pp.  9 - 15* 


i 
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For  two  consecutive  intervals  if  the  contributions  to  A(o>)  are  added, 


. - 


si^cntg)  for  the  first  interval  cancels  with 


"f(tlf| 

sin  (cut, 

CD  J 1 


) for  the 


second  interval.  Hence  the  sine  terms  drop  out  with  the  exception  of  the 

first  and  the  last  values  of  f(t).  Therefore,  if  t represents  the  initial 

value  of  t and  t the  last  value  of  t, 
n 


A(oj)  = 


f(t  J sin  (cut  ) 


^♦1 


Af  Acos((jot) 


Likewise, 


B(cd)  = 


f(t  ) - f(t)cos (ot  ) 


^X 


Af  Asin(ocit) 


The  previously  described  method  relies  on  the  assumption  that  the  function 
f(t)  is  linear  between  any  two  consecutive  points. 

Now  that  A(o>)  and  B(cjo)  are  known,  the  Fourier  transform  of  f(t) 


can  be  defined 


F(a>)  = n/a^cd)  + B2(cd) 


for  oo,  the  angular  frequency. 


The  Asin  cut  and  Acos  cut  were  computed  using  a recursive  method 
19 

suggested  by  R.  W.  Hamming.  y Since  the  time  interval  is  invariant, 

Asin(a>t  ) can  be  written  Asin(amAt),  and  Acos  (cat  ) can  be  written  Acos(omAt). 


^R.  y.  Hamming,  Numerical  Methods  for  Scientists  and  Engineers  (McGraw-Hill 
Book  Company,  Inc.,  1962),  pp.  72  - 73- 


The  first  step  in  the  recursion  is  to  define 


V = 0 
o 


> 


V 


1 


1 


> 


V = [2cos(cot)]v  . - V „ (m  = 2,5,  • ■ • ) 
m v m-1  m-2  v 7 


> 


and  then  show  that 


V 


m 


sin(mcot) 
sin  (cut) 


Since  the  above  equation  is  true  for  ra  = 0 and  1,  it  is  sufficient  to  show 


[2cos(cut)]v  , - V 0 
' m-1  m-2 


2cos(cot  )sin[  (m-l)cot] 
sin(cot) 


sin[  (m-2)  cut] 
sin  (cot ) 


sin(mcot)  + sin[  (m-2)cot]  - sin[  (m-2 )cot] 
sin(cot ) 


sin(mcot)  _ 
sin  (cot ) - vm 

Similarly  it  can  also  be  seen  that 

cos  (moot)  = [cos(t)]v  - V , 
v m ra-1 

One  further  consideration  when  computing  a Fourier  transform  is 
that  of  the  Nyquist  or  folding  frequency.  Given  a sampling  rate  At,  the 
sampling  theorem  requires 


2coAt  <1  , 

where  co  is  the  largest  angular  frequency  in  the  time  series.  As  will  be 
seen  in  the  following  section,  computation  of  the  Fourier  transform  for 
frequencies  near  co  has  an  undesirable  effect  on  the  transform. 


Input  for  FOURXFRM 
DATA  CARD 


DT,FREQ1,FREQEND,  DELTAF, IFPRINT  Y 4F10.5,I5  (FORMAT) 

INPUT  TAPE  FORMAT 

N 

Y(l)  Y(N) 

N 

Y(l)  Y(N) 

N 

Y(l)  •"  Y(N) 


END- OF- FILE 


program  FOURXFRm 

DIMENSION  Y (2b00) , C(300>,  DELTACOS(25UO) , 0£I_TASIN(?500) , 
l V (2500) * F ( 300) 

PI  = 3.141592654 
NMAX  = 0 
JF  s 0 
NL40  = 0 
V(l)  = 0.0 
V (2)  = 1.0 

READ  15.  D7,  FREQ).  FHEOEnD,  DELTAF,  IFPKINT  Y 
FRQ  = FREQEnD  ♦ 1.0E-08 
DO  100  I = 1.2500 

100  Y(I)  = 1.0 

101  IF  (JF  .EQ.  1)  108.  90 

90  BUFFER  IN  (1,1)  (N,N) 

N = N 

102  GO  TO  (102.103.108, 104)  , IJMTSTF(I) 

103  NMIN  = NMAX  ♦ 1 
NMAX  = NMIN  ♦ N - 1 

IF  (NMAX  .GE.  2500)  91,  105 

91  NMAX  = 2500 
JF  r 1 

GO  TO  105 

104  PRINT  10 
GO  TO  9999 

105  BUFFER  IN  (1,1)  (Y (NMIN) ,Y (NMAX) ) 

106  GO  TO  (106,101,107, 104) , UNITSTF(l) 

107  PRINT  11 
GO  TO  9999 

108  PRINT  12 
NOY  = NMAX 

DO  94  I s 1 , NOY 

IF  ( ABSF ( Y ( I ) ) .LT.  1.0E-40  .AND.  Y < I ) .NE.  0.0)  92,  94 

92  NL40  « NL40  ♦ 1 

IF  (NL40  .GE.  10)  93,  94 

93  IZERO  = I 
GO  TO  95 

94  CONTINUE 

95  DO  96  I * IZERO, NOY 
Y(I)  = 0.0 

96  CONTINUE 

IF  (IFPRINTy  .EQ.  1)  109,  110 

109  PRINT  13 

PRINT  14,  (Yd).  I = 1,NOY,500> 

PRINT  12 
PRINT  16 

110  J = 1 

FREQ  = FHEQl 
NM1  * NMAX  - 1 
U2  A « 0.0 
b * 0.0 

W s 2*  0*PI*FHtQ 
SINE  « SINF (**OT) 

COSINE  * COSF(W*DT> 

COSPROO  = (COSINE  - 1.0)* (COSINE  ♦ 0.5) 

ONEMCOS  * 1.0  - COSINE 
TOCOSM1  = 2. 0*COSINE  - 1.0 


I 
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CONS  * 1.0/ <»l*W*OT) 

M = 0 

A a A ♦ W*DT*SINF (W*Nm)*OT> *Y (NMAX) 

B a B ♦ W*OT*<Y(l)  - Y(NMAX)*COSF<rt*NMl*OT) ) 

113  M a M ♦ 1 

IF  (M  .EQ.  1)  114,  lib 

114  OELTACOS(M)  a -ONEMCOS 
OELTASIN(H)  a SINE 

GO  TO  118 

115  IF  (M  .EU.  2>  117,  116 

116  Y (M)  = 2.0*V(M-1)*COSINE  - V(M-2) 

117  DELTASIN(M)  a SINE* ( V (M) *TOCOSMl  - V<M-1)) 
DELTACOS  CM)  = 2, 0*V (M) *COSPROD  ♦ V (H-l ) *ONEMCOS 

118  A a A ♦ DELTACOS(M)*(Y(M*l)-Y(M) ) 

B a B * DELTASIN(M)* (Y (M,  1 )»Y (M) ) 

IF  (M.EQ.  NMl)  119,  113 

119  A a A*CONS 
B a 8*C0N5 

C(JJ  a SOHTF ( A*A  ♦ B*B) 

F(J)  a FREQ 

FREQ  = FREQ  ♦ DELTAF 

IF  (FREU  .LT.  FRO)  121,  122 

121  J a J ♦ 1 
GO  TO  112 

122  CMIN  = C < 1 » 

CMAX  = C(l) 

IMIN  = 1 

DO  128  I = 1»J 

IF  <C(I)  .GT.  CMAX)  125,  126 

125  CMAX  = C(I> 

GO  TO  128 

126  IF  (C(I>  .LT.  CMIN)  127,  128 

127  CMIN  = C<I> 

IMIN  = I 

128  CONTINUE 

DO  1029  I = 1»J 

cm  = cm*cm 

C ( I > = C'I)/(CMAX*CMAX) 

1029  CONTINUE 

DO  1030  I * 1»J 
PRINT  17,  F (I)  , C < I ) 

1030  CONTINUE 
CMAX  = 1.0 
CMIN  = C ( IMIN) 

PRINT  12 

CALL  MINMAX  (F(1),F(J),IFMIN,IFMAX) 

NOF  a IFMAX  - IFMIN 

CALL  MINMAX (CMIN,  CMAX,  ICMIN,  ICMAX) 

NOC  = ICMAX  - ICMIN 
IF  (NOC  .GT.  5)  129,  130 

129  YHEIGHT  = 1.0 
GO  TO  131 

130  YHEIGHT  a 2.0 

131  CALL  LAX(YHE1GHT,N0F,N0C) 

CONI  = IFMIn*2.0 

CON2  a 0.4 34 29* YHEIGHT 
CON3  * ICMIn*YHEIGHT 
LPEN  » 3 
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DO  132  K = 1»J 

XX  s 0.8685ft*ALOG <F (KJ ) - CONI 
YY  s CON2*AlOG(C(K) > - CON3 
CALL  PLOT (XX.YY.LPEN) 

132  LPEN  = 2 

CALL  PLOT (N0F*2. 0,0. 0,-3) 

10  FORMAT  ( 30H  PARITY  ERROR  ON  MAGNETIC  TAPE) 

11  FORMAT  (19H  EOF  IN  WRONG  PLACE) 

12  FORMAT  (1H1) 

13  FORMAT  (13X.2BHLISTING  OF  DATA  POINTS.  Y ( I ) ) 
1A  FORMAT  (6E20.10) 

15  FORMAT (AF 10*5* 15) 

16  FORMAT  (10X,AhFREO,lOX*I2HFOOH.  TRANS./) 

17  FORMAT  (E1A.A.3X.E19.5) 

9999  PRINT  12 

END 


SURROUTINE  MINMAX  (XMIN*  XMAX*  IMIN*  IMAX) 
I = 0 

1 IF  (XMIN  .LT.  10.0**(-I))  2*3 

2 1*1*1 
GO  TO  1 

3  IF  (I  .EQ.  0)  A.5 
A IF  (XMIN  .GT.  10. 0**1)  6.7 
6 1 = 1*1 
GO  TO  A 
7 IMIn  =1-1 
GO  TO  8 
5 IMIN  * -I 
81=  IMIN 

10  IF  (XMAX  .GT.  10. 0**1)  9.11 
9 1 = 1*1 

GO  TO  10 

11  IMAX  = I 

end 

SUBROUTINE  lax  ( YHEIGHT.NOF.NOC) 

CALL  PLOT (3. 0.-12. 0.-3) 

CALL  PLOT  (0.0. 0.5, -3) 

Y = 0.0 

DO  50  I = l.NOF 
A = 10.0**(I-1 ) 

DO  51  J = 1.10 
X = J 

X = 0 . 86858* ALOG ( X* A ) 

51  CALL  SYMBOL (X,Y,O.l»3»0.0,-2) 

50  CONTINUE 

CALL  PLOT  (0.0, 0.0, 3) 

X * 0.0 

CON  * 0.A3A29*YHEIGHT 
DO  52  I = l.NOC 
A « 10.0**<  T-l ) 

DO  53  J = 1,1) 

Y » J 

Y = CON* ALOG ( Y* A) 

53  CALL  SYMBOL (X.Y. 0.1.3, 0.0, -2) 

52  CONTINUE 
ENO 
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A 

B 

C(I) 

CMAX 

CMIN 

COSINE 

COSPROD 

DELTACOS(I) 

DELTAF 

DELTASIN(I) 

DT 

F(I) 

FREQ1 

FREQEND 

ICMAX 

ICMIN 

IFMAX 

IFMIN 

IFPRINT  Y 

LAX 

LPEN 

MINMAX 
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Glossary  of  Terms  in  FOURXFRM 

Accumulative  storage  location  for  sine  terms 

Accumulative  storage  location  for  cosine  terms 

Fourier  transform 

Maximum  Fourier  transform 

Minimum  Fourier  transform 

cos (o&t ) 

[cos(aAt)-l] • [cos(u£t)  + 0.5] 
cost  (N+l)afvt]  - costNdfit] 

Frequency  increment 
sint  (N+l)afit]  - sintNaAt] 

Time  increment 

Frequencies  for  which  Fourier  transform  is  computed 
Initial  frequency  for  which  Fourier  transform  is 
computed 

Last  frequency  for  which  Fourier  transform  is  computed 
Maximum  power  of  ten  of  Fourier  transform  range 
Minimum  power  of  ten  of  Fourier  transform!  range 
Maximum  power  of  ten  of  frequency  range 
Minimum  power  of  ten  of  frequency  range 
Parameter  to  determine  whether  or  not  to  print  input 
time  series 

Subroutine  to  plot  logarithmic  axes 
Parameter  used  in  plot  routine 

Subroutine  to  determine  maximum  and  minimum  powers 


of  ten  for  a sequence 
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NM1 

NMAX 

NMIN 

NOC 

NOF 

NOY 

ONEMCOS 

SINE 

TOCOSM1 

V(I) 

W(I) 

y(D 

YHEIGHT 


Number  of  points  of  input  time  series  minus  one 
Changing  parameter  determining  maximum  subscript  for 
each  block  of  buffered  data 

Changing  parameter  determining  minimum  subscript  for 
each  block  of  buffered  data 

Number  of  powers  of  ten  within  which  all  Fourier 
transforms  fall 

Number  of  powers  of  ten  within  which  all 
frequencies  fall 

Total  number  of  points  of  input  time  series 

1.0  - cos(o&t) 

sin(oAt) 

2.0cos(u&t)  - 1,0 

Parameter  used  in  recursion  for  DELTACOS(l) 
and  DELTASIN(l) 

Angular  frequencies 
Time  Series 

Parameter  used  in  plot  routine 


1 


L 
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IV.  EXAMPLE  RUNS  AND  RESULTS 

Four  examples  are  now  presented.  For  each  example  the  power 
spectral  density  function  V(ca)  and  the  time  interval  At  must  he  known.  The 
program  GAUSSIAN  is  then  used  to  obtain  an  approximation  to  the  impulse 
response.  FOURXFRM  is  next  used  to  obtain  the  Fourier  transform  of  the 

impulse  response.  The  square  of  the  transform  is  plotted  against  frequency 

20 

on  a log- log  plot.  The  examples  used  are  those  presented  by  Matuska. 

Because  of  the  manner  in  which  terms  of  the  w-sequence  are  used 
in  GAUSSIAN  to  form  vectors  for  the  solution  of  the  differential  equations 
(2.9)  and  (2.10),  the  impulse  time  series  was  formed  in  different  ways  for 
each  example.  If  the  degree  of  the  polynomial  Q(icn)  in  Eq.  (2.6)  was  n, 
then  the  impulse  of  magnitude  one  was  placed  at  w(n+l),  w(n+2),  • • • , 
w(2n)  and  the  impulse  response  computed  for  the  impulse  at  each  of  these 
times.  The  impulses  could  have  been  placed  at  w(kn+l),  w(kn+2),  • • • , 
w(kn+n),  k > 0,  with  the  same  resultant  impulse  responses.  It  is  seen  that 
n approximations  to  the  impulse  response  were  computed  for  each  of  the  four 
examples . 

The  theoretical  power  spectral  density  function  of  example  one  is 
given  to  be 

,T  / i Sto2  + 1 

V®)  = n 0 

cd  - 6o>  + 25 

which  has  zeros  at  the  points  (0,  1/3)  and  (0,  -1/3)  and  poles  at  the  points 
r,  r,  -r,  -r  for  r = (2,l).  Therefore  the  condition  0 < V^cn)  < <»  holds  for 


20, 


Matuska,  o£.  cit . , pp.  31  - 35- 
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all  real  03.  Also  V-^oa)  = V-^-cd)  and  V-^03)  -»  0 as  03  ->  ±°°  for  03  real.  With 
these  three  conditions  satisfied  V^(oo)  can  be  represented  in  the  form  of 
Eq.  (2.6) 

»<„),! 2 , 

I ( ico)  + 2(ico)  + 5 

p 

where  (ico)  + 2(io))  + 5 has  zeros  at  (-1,2)  and  (-1,-2).  Both  of  these 
zeros  lie  in  the  halfplane  with  Re(ios)  < 0. 


The  theoretical  power  spectral  density  function  of  example  two  is 


given  to  be 

k 

4 2 

„ / \ 03  - 70o>  + 1369 

vp(u>)  = -5 r 2 

03  + l4o3  + 4903  + 36 

It  can  be  shown  that  the  conditions  (2.5)  hold  such  that  ('03)  can  be 
represented 

P 2 

V (03)  = (iQ3)  +2(iy_)_+^7 

(103)5  + 6(io3)2  + ll(io))  + 6 

For  example  three  the  theoretical  power  spectral  density  function  Is 


„ / x 40003  + 1 

V,(03)  = -7T J7 0 

^ 03  l4o3  + 4903  + 3 6 


and  as  before  it  is  seen  that 


V?(03) 


I 20(iQ3)  + 1 

(i03)5  + 6(io3)2  + ll(io3)  + 6 


The  theoretical  power  spectral  density  function  of  example  four 


is  given  hy 


2 /.  2 2 > 
03  + (k  + c , 


V4(a>)  = 4Ak  -5 w ^ , 

us  + 2(k  -c  )ol>  + (k  +c  ) 

where  A = 2149.2,  k = O.658,  and  c = 2.71-  (cjd)  can  be  insured  to  have 
the  representation 

b (ioo)  + b 2 

V4(o>)  - 1 2 , 

a1(iaj)  + a2(icn)  + a^ 

where  b1  = 75-210999197,  = 53-521755647,  = 1.0,  a£  = 1.316,  and 

= 0.506405. 

For  each  of  these  four  examples  the  square  of  the  Fourier 
transform  of  the  impulse  response  was  found  to  agree  quite  well  with  the 
theoretical  power  spectral  density  function.  There  were,  however,  two  points 
at  which  discrepancies  did  occur.  One  of  the  points  was  at  the  extremely  low 
frequencies,  assuming  that  the  peak  frequency  did  not  occur  there.  For  those 
examples,  namely  two  and  four,  at  which  the  peak  frequency  was  the  lowest 
frequency,  there  was  no  discrepancy  between  the  theoretical  power  spectral 
density  function  and  the  square  of  the  Fourier  transform  of  the  impulse 
response.  In  examples  one  and  three  the  square  of  the  Fourier  transform  of 
the  impulse  response,  for  the  impulse  at  different  points,  bracketed  the 
power  spectral  density  function  at  these  low  frequencies. 

The  other  discrepancy  occurred  at  the  highest  frequencies.  As 
discussed  earlier,  the  frequency  range  is  bounded  above  by  (|  At)  where  At 


Ibid. , p.  34. 
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is  the  time  interval  between  points  of  the  input  time  series.  In  all 
examples  the  maximum  frequency  for  which  the  Fourier  transform  was  computed 
was  equal  to  j At.  This  fact  caused  a folding  back  of  the  Fourier  transform 
in  such  a manner  as  to  cause  the  square  of  the  Fourier  transform  of  the 
impulse  response  and  the  theoretical  power  spectral  density  function  to 
disagree  in  a small  manner. 

For  plotting  purposes  both  the  power  spectral  density  function 
and  the  square  of  the  Fourier  transform  of  the  impulse  response  were 
normalized  to  have  a maximum  value  of  one.  The  plots  of  these  four 
examples  are  shown  in  Figs.  1 through  4. 
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FREQUENCY  - Hz 
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APPENDIX 


The  white  noise  sequence  used  in  solving  the  differential 
equations  of  Franklin's  method  was  generated  by  taking  the  fractional  part  of 

9n  n = 1,  2,  3,  • • • 

where  0 is  a transcendental  number  greater  than  one.  The  fractional  part  of 
a number  $ is  represented  ($},  for  example  (3.14}  = 0.l4.  Franklin  has  shown 
that  a sequence  generated  in  this  manner  is  completely  equidistributed  for 
almost  all  0 > 1.^ 

Because  of  the  finite  word  length  of  a digital  computer  the 
"transcendental"  number  is  not  actually  transcendental.  Statistical  tests 
were  made  on  various  white  sequences,  and  in  all  cases  they  were  found  to  be 
equidistributed  with  a high  confidence  percentage. 

A question  of  interest  then  arises.  How  much  correlation  is  there 
between  sequences  generated  by  two  close  transcendental  numbers?  To  answer 
this  question  a program  was  written  to  obtain  the  cross-correlation  function 
between  (0n)  and  ((0+e)n),  for  n = 1,  2,  3,  * * * > N and  for  e arbitrarily 
small,  yet  still  within  the  machine  capability.  The  program  called  CORRELAT 
computes  the  cross -correlation  coefficients  for  ±10$  of  N lags,  where  N is 
the  total  number  of  points  in  the  white  sequence. 

No  theorem  was  proved  concerning  this  cross-correlation  analysis. 
Indeed  no  theorem  could  be  proved  because  of  the  fact  that  (0+e)  might  well 
be  zero  for  the  machine  capacity. 

22 

Joel  N.  Franklin,  "Deterministic  Simulation  of  Random  Processes," 

Math.  Comp.  1£,  (1963),  PP-  4 7 - 50. 
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An  example  run,  with  the  cross-correlation  function  shown  in 
Fig.  5,  was  chosen  as  follows: 

e = 2.71828183  2 e , 

-8 

e = 1.0  X 10  , 

N = 4266 

As  shown  in  Fig.  5,  there  was  no  significant  correlation  for  any 
of  the  lags  considered. 
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PROGRAM  COMPEL Al 

DIMENSION  X(4270).  Y (4270) .AM (402) ♦ C(402>»  IX(4),  IY(4),  M(400) 
I DEC  = 0 
ITEST  = 0 
CORRJ  = 0.0 
K = 0 
1STOP  = 0 
JSTOP  = 0 
XAVE  = 0.0 
YAVE  = 0.0 
DO  260  J * 5»b 
NEOF  * 0 
READ  10,  1EOF 
IF  ( IEOF  . E(J.  0)  260.  80 
BO  I IEOF  a 2* IEOF 
90  ITEST  a o 

BUFFER  IN  'J,))  (IHUF.IMUF) 

100  60  TO  (100.lf0.130.110)  ♦ UMTSTF(J) 

110  PRINT  12 
GO  TO  9999 
130  NEOF  = NEOF  ♦ 1 

IF  (NEOF  .EQ.  1 1 E OF ) 260,  90 

140  IBUF  = I BUF 

IF  (IBUF  .LT.  500)  141.  142 

141  ITEST  a 1 

142  BUFFER  IN  (J.l)  ( X ( 1 ) * X ( I JUF ) ) 

150  60  TO  ( 150.160. 160.110)  . IJNITSTF(J) 

160  IF  (ITEST  .EQ.  0)  90*  170 
170  BUFFER  IN  (J.l)  (IOEC.IOEC) 

1 BO  GO  TO  (180.200,190.110).  UNllSTF(J) 

190  NEOF  = NEOF  ♦ 1 
GO  TO  170 

200  BUFFER  IN  (j.l)  (L.L) 

210  GO  TO  (210.  230.  230,  110),  UNITSTMJ) 

230  BUFFER  IN  (j.l)  (M ( 1 ) ,M (400) ) 

232  GO  TO  (232.90,90,110),  UNITSTF(J) 

260  CONTINUE 

265  BUFFER  IN  (6,1)  (I, I) 

270  GO  TO  (270.  2B0»  280,  110),  UNITSTF(6) 

280  IF  (I  .EM.  500 ) 300.  290 
290  ISTOP  = 1 
300  K1  = 500*K  ♦ 1 
I a I 

K I a 500*K  ♦ 1 

BUFFER  IN  (6,1)  ( X (M  ) , X (K  I ) ) 

310  GO  TO  (310,  320.  320,110).  UNITSTF<6) 

320  BUFFER  IN  (5,1)  (J.J) 

330  GO  TO  (330,  340,  340,  110),  UNITSTF(5> 

340  IF  (J  .EO.  500)  3b0.  35o 
350  JSTOP  a 1 
360  KK1  a 500*K  ♦ 1 
J a J 

KKI  a 500*K  ♦ J 

BUFFER  IN  (5*1)  ( Y ( K.K  l).Y(KKl)) 

370  GO  TO  (370.380.380, 1 10)  « UMTSTF(5) 
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360 

390 

400 

410 

420 


430 


440 

450 

460 

470 

460 

490 

500 


10 

11 

12 

13 

14 

15 


K = K ♦ 1 

IF  ( I STOP  .EG.  ft  .AND.  JSTOP  .£0.  </)  ?e>5*390 

IF  (I-J)  400*  400*  4 1 U 

N r I ♦ (K-l)*500 

60  TO  420 

N = J ♦ (K-l ) *500 

DO  430  K = 1*im 

XAVE  = XAVE  ♦ A(K> 

YAVE  = YAVE  ♦ Y ( l\ ) 

CONTINUE 
XAVE  = XAVE/N 
YAVE  * YAVE/N 
PRINT  17*  X AVt , YAVE 
DO  440  K = 1*N 
X (K)  = X ( K ) - XAVE 
Y (K  ) = Y ( K ) - Y AVt 
CONTINUE 

READ  11*  ML,  INC 

MM  s 2* I AbS (ML ) / INC  ♦ 1 

DO  490  L = 1**M 

ML  = ML  ♦ INC 

IF  (ML)  450  » 450*  460 

NN  = N ♦ ML 

K = 1 

60  TO  470 

NN  s N 

K = ML  ♦ 1 

DO  460  I = K * NN 

J * 1 - ML 

CORRJ  = CORRJ  ♦ X ( 1 ) * Y ( J ) 

CONTINUE 

C(L)  = CORRJ/ (N-IAHS(ML) ) 

CORRJ  = 0.0 

CONTINUE 

AM ( J ) = -200.0 

DO  500  I = 2 * MM 

AM  ( I ) = AM(T-l)  ♦ 1.(1 

PRINT  13,  MM 

PRINT  14*  (C ( J) * J = 1 ,MM> 

IX(1)  = 4HLAb 
IY(1)  = 4HC 

CALL  PLOT  <3. 0,-12. o*-3) 

CALL  PLOT  (0.0*0. 5,-3) 

MM1  = mm  ♦ 1 
MM2  = MM  ♦ ? 

c(mmd  = -o.ioo 

C (MM2)  = 0.025 
AM(MMl)  = -200.0 
AM(MM2)  = 20.0 

CALL  AXIS  (0. 0*0. 0*IX«-4«?0. 0*0. 0,-200. 0*20. 0,0. 1 *1. 0*1  ) 
CALL  AXIS  (0. 0*0. 0»IY, 4, 6. 0*90. 0,-0.10,0.025, 0.1, 1.0,1) 
CALL  LINE  <AM,C,MM,1, o» 0* 0. OB, 1 ♦ 0 . 0) 

FORMAT  (12) 

FORMAT  (1615) 

FORMAT  ( 1H1 ,2 1H  PARITY  ERROR  ON  TAPE) 

FORMAT  (35X16M  CORRELATION  FOR,I4,5H  LAbS//) 

FORMAT  (6(4XF 12.H) ) 

FORMAT  (2 (4Xt 1 7.10) ) 
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I 
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16  FORMAT  ( 1 Ml ) 

17  FORMAT  (13H  X AVtXAOt  = »F17.10»13H  Y A VfcwAGl  = iHM'  ) 
9999  PRINT  16 

ENO 


I 

I 
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C(I) 

CORRJ 

IBUF 

IDEC 

IEOF 

INC 

ITEST 

ISTOP 

IX(I) 

IY(I) 

JSTOP 

ML 

NEOF 

X(I) 

XAVE 

Y(I) 


Cross-correlation  coefficient  for  lag  I 

Temporary  storage  for  correlation  coefficients 

Number  of  points  in  each  block  of  input  tape 

Number  of  points  in  last  block  of  input  tape 

Number  of  sequences  on  tape  prior  to  sequence  of  interest 

Number  of  points  skipped  between  each  lag 

Parameter  to  determine  location  of  data  on  input  tape 

Parameter  to  determine  whether  all  of  the  X(l)  sequence 

has  been  buffered  from  tape 

Abscissa  label  for  plots 

Ordinate  label  for  plots 

Parameter  to  determine  whether  all  of  the  Y(l)  sequence 
has  been  buffered  from  tape 
Number  of  lags 

Number  of  end-of-file 's  on  tape  prior  to  sequence  of 

interest 

Input  sequence 

Average  of  X(l)  for  lag  I 


i 
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YAVE 


Input  sequence 

Average  of  Y(l)  for  lag  I 
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